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Abstract. We study the statistical mechanics of a model describing the coevolution of 
species interacting in a random way. We find that at high competition replica symmetry is 
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and we compare our findings with accurate numerical simulations. 
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1. Introduction 

The replicant models study the coevolution of sets of interacting species able to re- 
produce themselves: they have a huge number of applications in biologic and optimization 
problems [1-4]. In this paper we study a non-deterministic evolution: we consider a 
system of replicants which evolve with random interactions. 

The model, introduced by S. Dicdcrich and M. Opper in [5], is defined as follows. 
Given species, let Xi/N be the concentration of the i-th family in the system. The real 
variables {xi G IR, i = 1, . . . , A^} are then subject to the constraints 

N 

J2^i = N, Xi>0 yi = l,...,N. (1.1) 

1=1 

The interactions between different species are described through a fitness functional 
Fj{x) that must be maximized at equilibrium. Typically, Fj is chosen as a quadratic 
function of the concentrations, that is equivalent to take into account only pair interactions 
between the species: 

N N 

Fj{x) ^ -nj[x] ^ - ^ Jij Xi Xj - a^xf, (1.2) 

i<j=l i=l 

where the parameters {Jij} are chosen at random from the Gaussian probability distribu- 
tion 

m.) = yj^exp(^-^j (1-3) 

like in the Sherrington-Kirkpatrick model of spin glasses [6, 7]. The control parameter a 
has the aim of limiting the growth and the supremacy of one single species: for big values 
of a, the growth of all the species is strongly limited by the factor axf; in that case, the 
random interactions become negligible and the equilibrium configuration is 

xf-l Vi = l,...,7V, (a>J) (1.4) 

almost independently from the interactions between the species. Instead, for small values 
of a, the pair interactions play a central role and a few species prevail among the others. 
Analytically, this model differs from the SK spin glass in that we impose the constraint 
(1.1): the spins are then allowed to take any real value, but the total magnetization is 
fixed. 

In section 2 we show how it is possible to solve the random replicant model within 
the replica formalism. In sections 3 and 4, we analyze the replica symmetric solution and 
its stability and in section 5 we perform the first step of the hierarchical replica symmetry 
breaking. The biological applications of the results are found in the limit T — > 0+ because 
the fitness functional Fj introduced in the last section is, a minus sign apart, the low 
temperature limit of the free energy. 
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The study of the stabihty of the rephca symmetric solution will show that, at zero 
temperature, the replicant model exhibits a phase transition to a glassy phase when a 
crosses a certain value ac- The replica symmetry breaking which occurs in the glassy 
phase (a < ac) implies the breakdown of the ergodicity of the system: when a becomes 
small, the evolution of the system depends strongly on the initial conditions, and in general 
we will not be able to make any precise prediction on the equilibrium state of the system. 

From the biological point of view, the glassy phase is the unstable one: in the high 
a phase, a single equilibrium state exists, and the system is able to recover its equilib- 
rium configuration after any external change of the concentrations of its elements; on the 
contrary, in the glassy phase, the same perturbation can change drastically the final con- 
figuration of the system, if it is led to a different ergodic region of the phase space. Here 
however we study only the properties of the statics associated to Hamiltonian (1.2) and 
we do not consider the dynamics of a system leading to this equilibrium distribution. 

2. The Random Replicant Model: analytical solution 

Now we derive the expression for the quenched free energy density of the random 
replicant model. In this and the next section, we follow closely [5]. The evolution of the 
system is ruled by the Hamiltonian (1.2); averaging over all the possible choices of the 
{Jij}, the quenched free energy of the system is given by 

-pNf = I JJdJi,- P(Jy) InJ^exp {-mj[x\) • (2.1) 

To compute (2.1), we use the replica method [7 - 12] introducing a set {Aq, a = 1, . . . , n} 
of Lagrange multipliers which ensure the normalization condition (1.1) in each of the n 
replicas. With standard calculations [13], we arrive at the following expression for /: 



—f3f = lim max 
n^0+ Q,\ 



■ - + - A„ + - In IVn exp L(Q, A, x) 



n ^ — ' ' n ^ — ' n 



(2.2) 



where 

L(Q, A, x) := -/3a XI ^« ~ XI '^"^^ + X QajXaX^ ; (2.3) 

a a ay 

Q and A are respectively an n xn matrix and an n-dimensional vector; {xa, ct = 1, . . . , n} 
is a new set of real positive variables, and Tr„ denotes the integral over all possible values 
of the XaS. 

From (2.2) and (2.3) we find that the stationarity equations for / are: 



PJ TrnXaXjexpL{Q,\,x) 

Qay = ^ T rr^ \ — ^ — Vo;, 7 = 1, . . . , n; 

2 TVn expL(Q, A, cc) 

Tr„.c<, cxpL(Q, A, a;) 

1 = — — ^ ^ — r — Vq! = 1, . . . , n. 

Tr„ expL(Q, A, x) 



(2.4) 



The remaining sections are devoted to the study of the solutions of these stationarity 
conditions. 
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3. Replica symmetric solution 



Both the free energy density and the stationarity equations above are invariant under 
the action of the group 5„ of permutations between the n rephcas. This imphes that at 
least one of the solutions of (2.4) is invariant under Sn, so that the first ansatz that is to 
be tried is certainly the symmetric one, which is given by: 



Qcxf — Q ^a-y ~l~ t 



and 



Aa — A. 



(3.1) 



Introducing (3.1) into (2.2), and denoting by /rs the resulting the free energy density we 
have, after the manipulations described in [13], 



-/5/rs = max 

q,i,\ 



g2 + 2(3 J qt-(3J X- In 



dxexp £-Eis{q,t, X,x,z) 



where 



C^s{q, t, A, X, z) := —j3J (a — g)x^ — {2z\/ 1 — X)x 
t := t/{pj), X := X/{PJ), a := a/ J, and we have introduced the notation 

1 f + °^ 2 

G{z) := / d^e-^ G{z). 
In a similar way, the stationarity equations become: 



{x)l = 2t, 



2q_ 
PJ'' 



where 



( G{x) ). 



+ 00 



dx G{x) exp Cjis{q,t, X,x, z) 



r+oo 



dx expjC^_siq,i,X,x,z) 



(3.2) 



(3.3) 



(3.4) 



(3.5) 



(3.6) 



The low temperature limit of the symmetric solution, first studied by Diederich and 
Opper [5], is particularly interesting because it allows us to prove analytically the exis- 
tence of a second order transition to a glassy phase, as we will show in the next section. 
Introducing the parameter 

■ ^ (3.7) 



z := 



2Vt ' 
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the stationarity equations (3.5) become: 



4:q{a-q) 




dz e 



< 4 t (2g - a) 



A; 



(3.8) 



V 



2(a - q) 




dz e 



leading to /rs = 2Jt(a — 2q). Figure 1 shows how g, i, A, and /rs behave as functions 
of a in this limit. We give also the approximate expressions of these parameters in two 
particularly interesting cases: the "classical" regime (a ^ 1), and the neighbourhood of 
the critical point Sc = 1/ -\/2. 



Figure 1. Numerical solutions of the replica symmetric equations in the low temperature 
limit 
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In the former region the equihbrium configurations become trivial, with x^^ ~ 1, Vz. 
The rephca symmetric solution, which we will prove to be stable in this region, predicts: 



A = -a-A/a2-l + C>(e-"'), 
and the free energy density becomes 

/rs = \ (a+ ^/^^) +0(e-«'). (3.10) 

Instead, the latter is the transition point to the glassy phase, as we will show below; in its 
neighbourhood we have: 

9 = ^ 2~*^^ ~ '^^^^ ~ ^^^^ + ^^^^ ~ 

i= I -A/27r(7r-2)(a-ac)+7r2(37r-8)(a-ac)^ + o((a-ac)^), ^^^^ 
A = -27r(a - ttc) + 2V2n{n - 2)(a - ttc)^ + o((a - acf), 
/rs = ^(^ - ^c) - V2n{TT - 2) (a - Oc)^ + o((a - Oc)^). 

Finally, figure 2 shows the numerical results that we obtained for the order parameters 
q and t by solving equations (3.5) for different finite values of /3. 



Figure 2. Numerical solutions of the replica symmetric equations at finite temperature 
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4. Instability of the symmetric solution 

In the preceding section we have shown that equations (2.4) admit a symmetric so- 
lution, but we must also check the Hessian of the free energy to determine whether our 
solution is a minimum of / or just a saddle-point. To find the eigenvalues of the Hessian 
we generalize the calculus made by De Almeida and Thouless [9] for the SK model of spin 
glass: let 

Qaf = iq + Sqa) Saj +t + Staf, with | ^f"" _ ^' . 

ULq,j Uijce, (4.1) 

If we denote by 6$, the vector {d\; Sq; St) and we substitute (4.1) in (2.2) we obtain, 
after some tedious calculations [13], that the second order term in the expansion of / in 
terms of is given by —(352f = \ 5^ ■ M. S^, where Al is a real symmetric matrix with 
the following fourteen different types of elements: 

A := Msx^sx^ = (WVz - MT^) , 

B:= Msx^sx,= (Wz-J^'), 

C := Msx^s,^ = -/3J {Wh - Wil , 

D := Msx^s,, = -PJ {WWu - Wu ¥u) , 

E := Msg^s,^ = -2 + p^J^ ((^ - (^') , 

F := Ms.^s,, = -p'j' (wn - wu') , 

) , ^ (4.2) 

H := Msx^5t-,s = -pj [i^n - i^n (^).) , 

K := Mst^,5u, = -2 + P^J-" (F)! - W^) , 
K' := Mst^^st,^ = P'J' - W') =K + 2, 

L := Mst^^st,, = p'j' (yhW. - W') , 

M := Mst^^st,, = P'J' (W - M!') • 

Furthermore, A4 has three different types of eigenvectors: 
(i) symmetric eigenvectors of the type 

S^={£,...,£;p,...,p;T,...,Ty, (4.3) 
6 



P Biscari & G Parisi 



(ii) 1-asymmetry eigenvectors, with: 

5Xa = 



£i if a = a 
£o otherwise, 



(iii) 2-asymmetries eigenvectors, of the type 

5Xa = 



r. I Pi if Q! = (5 ^ A A\ 

Sqa = V ■ (4-4) 

' po otherwise, ^ ^ 

Ti if q; = (5 or 7 = q; 

To otherwise; 



£i if a = CK or a = 7 
£o otherwise. 



r _ j Pi ifa = Q:orQ! = 7 

~ \ Po otherwise, (4-5) 

{T2 if q;7 = 0:7 or 0:7 = 70: 
To if q; 7^ (5, a 7, 7 7^ Q! and 7 7^ 7 
Ti otherwise. 

The eigenvalues of M. must be negative in order to ensure the stabihty of the symmet- 
ric ansatz. The biggest among the eigenvalues associated with the families above comes 
from the 2-asimmetries eigenvectors, and is given by [13] 



/xer = -{K + K'-2L + M) = -l + f3\{x^), - {x)lf. (4.6) 
In the low temperature limit //cr can be easily computed, and is equal to 

l^cr = • (4.7) 

a — q 

In particular, as figure 3 shows, ficr becomes positive when a < a^. 

A*cr = 7r(a — he) — -\/27r(7r — 2) (a — hcf' + o((a — Sc)^). (4-8) 

5. Replica symmetry breaking 

Having proved the instability of the symmetric solution, we must now search a more 
general ansatz to describe the system when a < a^. To obtain it, we will follow the 
guidelines of the hierarchical ansatz of spin glasses [10, 11, 12]. In this paper we study 
only the first step of the replica symmetry breaking, testing order parameter matrices of 
the type 



\ if Int (^j ^Int 



Qaj — * 



a 



t + r if Int y-j = Int ) but a 7^ 7, 
. q + t + r ifQ: = 7. 



1\ ^ (5-1) 
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Figure 3. Critical eigenvalue at zero temperature. The symmetric solution becomes un- 
stable when it is positive. 



This ansatz can be improved by iterating the breaking scheme in all the blocks introduced in 
the first step, but we will see that even a single breaking improves drastically the symmetric 
predictions. We recall that, in the limit n — > 0+, the hierarchical parametrization can be 
written in terms of an order parameter function Q{x), defined in the interval x G [0, 1], 
which, at this point of symmetry breaking, is equal to 

In (5.2) we have omitted the diagonal term containing q (corresponding to (5(1)), because 
it involves the term in the Hamiltonian that contains a and it can always be treated 
separately. The introduction of the breaking parameters rj and r changes the free energy 
density as follows: 



-/?/h = max 

q,t,r,X,7] 



- {q + t + rf - {r]-l){t + rf + ijr + A+ 

(5.3) 



H — / d^ — In / d^r — 1^ \ / dx exp>CH(Q, ^, A, x, 2, ^r) I 
with 

>Ch(9, r, A, X, z, Zr) := -/?J(a - q)x'^ + {2z^/pji + 2zr\f^Jr - A)a;. (5.4) 
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The stationarity equations related to /h become: 



(x) 



(3J 

PJJ 
2 V 



{x) 



{z,Zr) 



(x) 



{z,Zr) 



{x) 



{z,Zr) 



(5.5) 



Q = 



2 V 



{z,Zr.) 



(x) 



{z,Zr) 



ri^r (r + 2t) = 





p + OO 






In / da; >Ch 


— In / dZr — ■i=^V{z^Zr) 1 




. Jo 


z J-oo VTl" / 



where 



+00 \ n 

da; >Ch(9, A, a;, 2;, Zr) 



{z, Zr) := (^J^ 

r+oo 

/ da; • jCii(q,t,r,X,x,z,Zr) 
\ Jo 



{z,Zr) 



r+oo 

/ da; jCh (q, t, r, A, X, z, Zr) 
Jo 

r+oo ^-zl 

/ dZr — j=- ■ Viz^Zr) 
J —00 



(5.6) 



/. 



+00 



dZr 



7T 



P{z,Zr) 



Solving numerically equations (5.5) in the low temperature limit, we find that the product 
rjr of the two breaking parameters remains finite in the /? — > +00 limit, and that it becomes 
different from zero as soon as a < Sc, as it is shown in figure 4. 

In figure 5 wc show one of the results that we found in the numerical simulations that 
we have performed on this model, and that we will describe in more detail elsewhere: the 
triangles represent the free energy density obtained from the simulations at zero tempera- 
ture, the continuous line corresponds to the replica symmetric prediction, and the dashed 
line illustrates the broken symmetry results. Figure 5 clearly shows how the first step of 
the replica symmetry breaking improves the symmetric predictions, even if it fails when a 
goes to zero. 

To conclude the study of the replica symmetry broken solution we will now show that, 
in the low temperature limit, ij goes to zero as C(/?~^), so that the breaking parameter r 
scales as t and A do, i.e. that r = f3r with f finite as T goes to zero. We will prove this 
result near the critical value of a where we have a better analytic control. To this end, we 
push our expansion of / to the third order in 5\ 6q, 6t, obtaining 



(2) 



-/3/ = -/3/rs+ lim 

n->0+ 



+ /(3) 



n 



(5.7) 
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Figure 4. Numerical solutions obtained for the product of the breaking parameters r] and 
r at zero temperature. 



Figure 5. Improvement led by the 1-step replica symmetry broken solution in the prediction 
for the free energy at zero temperature. 

The second order term /^^^ was studied in the preceding section; considering the third 
order term as a functional of the order parameter function Q{x), and neglecting higher 



10 



P Biscari & G Parisi 



order terms in n, the stationarity equation 

S 



5Q{x) 



f[q] = 0, 



(5.8) 



that must be verified Va; e [0, 1], can be given the form [13]: 



+ 



+ 2i(3^{x^), - 3{x^)Ax)z + 2{x)l) {{x^), - {x)l) + 



- 2pl3'{{x^), - 2(x3),(x), - + 2(x)2(,x2),) ((,t2), - {x)l) + 



+ 2p^Q{x)i2{{x^), - {x)iyx - ((a;3), - {x^),{x),){{x^), - 5{x^),{x),+ 



(5.9) 



+ A{x)l)+A{x)l{{x^),-{x)iyj + 

+ 4/?3 dx' g(x')(F>^ - s{x^)z{x), + 2Wz) {WVz -Wz)+ 

Jo 

+ 4(3' [ dx'g(a;')(F)I-M!)'}=0, 



where 



Q{x) 



dQ 
dx 



(5.10) 



This expression can be derived again with respect to x to obtain a necessary condition 
for the equihbrium: 



Q{x) = or 
-4/33 



2x1 - {x)l - {x'), - 3{x^).{x), + 2{x)l 



(5.11) 



This is exactly what we were fooking for: the order parameter function Q{x) must be 
a constant in all x e [0, 1], except for 



X = T] 



{x'),-3{x^),{x), + 2{x)l 



2 {x^).-{x)] 



(5.12) 



where a jump can happen. Note that in Ising spin glasses with two spins interaction no 
solution of this type can be found, while a similar phenomenon happens in the Potts model 
with p components, when p > 4. 
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The integrals in (5.12) can be computed in the low temperature limit, leading to [13]: 



where z is defined in (3.7), and the numerical constant cn can be easily evaluated, and is 
equal to 0.0167066... Equation (5.13) shows that the scaling behaviuor of r] is precisely 
ry = 0{P~^) when (3 +oo. 

6. Conclusions 

We have shown that in the replicant model replica symmetry is broken. The pre- 
dictions based on one step replica symmetry breaking are in better agrement with the 
numerical data than those coming from exact symmetry. Contrary to what happens in 
most of the cases also the one step replica symmetry breaking is not able to capture the 
behaviour of the system in the limit of very small a. It would be rather interesting to 
obtain the results from full replica symmetry breaking in this region. This task should not 
be impossible using the techniques of [13]. 
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